Childhood exposures to environmental chemicals and neurodevelopmental outcomes in congenital heart disease

Background Children with congenital heart defects have an increased risk of neurodevelopmental disability. The impact of environmental chemical exposures during daily life on neurodevelopmental outcomes in toddlers with congenital heart defects is unknown. Methods This prospective study investigated the impacts of early childhood exposure to mixtures of environmental chemicals on neurodevelopmental outcomes after cardiac surgery. Outcomes were assessed at 18 months of age using The Bayley Scales of Infant and Toddler Development-III. Urinary concentrations of exposure biomarkers of pesticides, phenols, parabens, and phthalates, and blood levels of lead, mercury, and nicotine were measured at the same time point. Bayesian profile regression and weighted quantile sum regression were utilized to assess associations between mixtures of biomarkers and neurodevelopmental scores. Results One-hundred and forty infants were enrolled, and 110 (79%) returned at 18 months of age. Six biomarker exposure clusters were identified from the Bayesian profile regression analysis; and the pattern was driven by 15 of the 30 biomarkers, most notably 13 phthalate biomarkers. Children in the highest exposure cluster had significantly lower adjusted language scores by -9.41 points (95%CI: -17.2, -1.7) and adjusted motor scores by -4.9 points (-9.5, -0.4) compared to the lowest exposure. Weighted quantile sum regression modeling for the overall exposure-response relationship showed a significantly lower adjusted motor score (β = -2.8 points [2.5th and 97.5th percentile: -6.0, -0.6]). The weighted quantile sum regression index weights for several phthalates, one paraben, and one phenol suggest their relevance for poorer neurodevelopmental outcomes. Conclusions Like other children, infants with congenital heart defects are exposed to complex mixtures of environmental chemicals in daily life. Higher exposure biomarker concentrations were associated with significantly worse performance for language and motor skills in this population.

Introduction Improved survival after surgical repair of congenital heart defects (CHD) has led to recognition that neurobehavioral disability, including mild cognitive dysfunction, impaired motor skills, language difficulties, and impaired attention and executive function, is the most common and potentially most disabling long-term adverse outcome of infants with CHD [1]. Despite the efforts of multiple investigators and clinicians, there has been minimal progress in improving early neurodevelopmental outcomes for these children over the past 20 years [2]. For many years, studies of neurodevelopmental outcomes for children with CHD were focused on cardiac surgery; specifically on intraoperative management strategies [2]. This focus was understandable inasmuch as these strategies are potentially modifiable and thus there were potential opportunities to improve outcomes. However, recent studies have shown that innate patient factors (prematurity, brain dysmaturity, genetic anomalies, maternal education, etc.) eclipse intraoperative factors as predictors of worse neurodevelopmental outcomes [1,2]. Moreover, multiple studies have demonstrated that currently recognized risk factors explain only a small portion of total variation in neurodevelopmental outcomes [2]. A better understanding of other modifiable, yet under recognized factors, is critical in not only identifying risk profiles, but also in guiding development of more individualized therapeutic strategies.
There is growing evidence of the adverse effects of even low level, early life exposures to environmental chemicals on neurodevelopment (e.g., metals, tobacco smoke, pesticides, industrial chemicals [phthalates, phenols, parabens]) [3,4]. Many of these chemicals are common in the environment and are known to be endocrine disrupting compounds (EDCs) [3,[5][6][7][8][9][10]. Household dust can contain lead (Pb), phthalates, and pesticides [11,12]. Phthalates, phenols, and parabens are present in many common consumer and personal care products [13]. In some cases (e.g. Pb), there is no known "safe" threshold for exposure to these environmental chemicals, with infants having increased exposure vulnerability and biologic susceptibility compared to other sub-populations. The early postnatal period through the first two years of life is a particularly important time for brain development [1,4]. The infant brain is in a state of rapid growth; thus, impairment of brain cognitive function may arise from relatively low levels of exposure [3,4]. Many previous epidemiologic studies have focused on exposure to one neurotoxicant at a time [4]. However, environmental chemical exposures occur as complex mixtures of chemicals from multiple chemical classes [4]. Several recent studies have suggested that exposure to environmental chemical mixtures may have cumulative or possibly synergistic (or antagonistic) adverse effects on neurodevelopment [4,[14][15][16][17].
We have previously shown that, in the perioperative period, neonates undergoing cardiac surgery with cardiopulmonary bypass (CPB) are exposed to high-levels of phthalates and phenols from medical devices [18]. Cyclohexanone is a major industrial solvent used in the production of medical devices [19]. Everett and colleagues found substantial exposure to cyclohexanone in neonates undergoing CPB [19]. Higher exposures were associated with worse performance for cognitive and language skills at 12 months of age. However, to our knowledge, no previous study has investigated the impact of exposure to environmental chemicals (in isolation or as part of a mixture) encountered during daily life on neurodevelopmental outcomes in the CHD population. It is important to realize that in the infant CHD population, these chemical exposures co-occur with other neurodevelopmental risks (e.g., brain dysmaturity, hypoxia, cardiac surgery, etc.), which may exacerbate the risk of neurotoxicity [1,4]. In this context, children with CHD are a particularly at-risk and vulnerable population.
The purpose of this study is to determine if early life chemical exposures are potentially important factors which are not currently included in our understanding of the mechanisms underlying neurodevelopmental disability in children with CHD. We assessed a large number (n = 30) of exposure biomarkers of neurotoxicants and EDCs, including Pb, mercury, environmental tobacco smoke (cotinine, nicotine), pesticides, phthalates, phenols, and parabens.

Patients and methods
We conducted a prospective, observational cohort study investigating the effect of childhood exposure to mixtures of environmental chemicals on neurodevelopmental outcomes at 18 months of age after cardiac operations in newborns. Inclusion criteria were: (1) infants with CHD and an expected operation with cardiopulmonary bypass (CPB) at age 44 weeks postconception or younger. Exclusion criteria were (1) presence of an identified genetic syndrome, (2) major extracardiac anomaly, or (3) language other than English spoken in the home. The Children's Hospital of Philadelphia (CHOP) Institutional Review Board and the University of Florida Institutional Review Board approved the study. Written informed consent was obtained from the parent or guardian. The involvement of the Centers for Disease Control and Prevention (CDC) laboratory did not constitute engagement in human subjects research.
Patients were enrolled between September 1, 2011 and August 31, 2015. Operations were performed by 1 of 4 cardiac surgeons with a dedicated team of cardiac anesthesiologists. Deep hypothermic circulatory arrest (DHCA) was used at the surgeon's discretion. Modified ultrafiltration was performed in all patients. Patient-related variables (e.g., age at testing, sex, race, anthropometric measures) and peri-operative variables (e.g., age at surgery, bypass support times, hematocrit, length of stay) were collected from patient records. Newborn length, weight, and head circumference were measured, and the z-score for each growth measurement was calculated using World Health Organization (WHO) standards for full-term infants and the Fenton growth chart for preterm infants [20,21]. Body mass index (BMI) and z-scores were calculated using WHO standards for full-term infants and the Olsen growth chart for preterm infants [22]. Maternal education, socioeconomic-related variables that comprise the Hollingshead socioeconomic status (SES), and ethnicity were determined through parental report [23]. A blood sample was obtained from each patient for genetic analysis. Whole exome sequencing was performed on subjects and all consented parents. Details are included in the Supplemental Files. Variant calls were queried for single nucleotide polymorphisms to determine apo-lipoprotein E (APOE) genotype according to the following classification of bases at rs429358 and rs7412, respectively on chromosome 19: C,T, ε1; T,T, ε2; T,C, ε3; C,C, ε4.
A study follow-up visit was conducted at 18 (mean) ± 1 (SD) months of age. Growth parameters (weight, height, and head circumference, BMI) were measured, and z-scores derived using WHO standards. Collection of spot infant urine samples was performed by placing a cotton ball in the diaper. The cotton balls were placed into a 3 to 5 mL polypropylene syringe (Becton-Dickinson, Franklin Lakes, NJ), and the urine was transferred into a polypropylene Cryovial (Simport, Beloeil, QC, Canada). Specific gravity (SG) was measured at room temperature using a handheld refractometer. Urine and negative control samples were frozen at -20˚C and shipped overnight on dry ice to the CDC. At the CDC, urinary biomarkers were quantified in infant urine using solid phase extraction coupled to high performance liquid chromatography-isotope dilution tandem mass spectrometry after enzymatic deconjugation to hydrolyze urinary conjugates and following standard quality assurance/quality control procedures as described in detail previously [24][25][26][27]. The target chemical biomarkers measured and their acronyms are shown in Table 1. The limits of detection (LOD) were 0.1-1.7 μg/L, depending on the biomarker. A venous blood sample was obtained for measurement of lead, mercury, nicotine, and cotinine. Lead levels for Pennsylvania residents were performed by the CHOP Clinical Laboratory. Measurements for out of state subjects were performed by ARUP Laboratories (Salt Lake City, UT). Total mercury measurements for all subjects were performed by ARUP Laboratories. Nicotine and cotinine were measured by ARUP Laboratories and NMS Labs (Horsham, PA).
Patients were evaluated by a genetic dysmorphologist. Neonatal recognition of dysmorphic features may be difficult; therefore, some patients for whom the diagnosis of a genetic syndrome was made at a later evaluation were enrolled. Genetic testing was performed as clinically indicated. Results of the genetic evaluations were classified as normal if no genetic or chromosome abnormality was demonstrated, abnormal if a specific diagnosis was confirmed, and suspect if there was evidence of a genetic syndrome that could not be confirmed. Neurodevelopmental outcomes were assessed using The Bayley Scales of Infant and Toddler Development-III; which provide composite scores for cognitive, language and motor skills (μ = 100 with σ = 15) [28]. Higher scores indicate better skills.

Statistical analysis
Data analysis proceeded in three phases: a descriptive analysis phase to fully characterize the study population; a modeling phase in which selected neurocognitive measures were regressed onto a series of established risk factors to identify influential covariates for subsequent modeling; and, finally, a health effects modeling approach which explicitly evaluated the effects of early-childhood chemical mixtures on neurodevelopmental attainment. All data were analyzed in SAS (v9.4) or R (v 4.1.2).

Descriptive statistics
Measures of central tendency, variability, and univariate associations with neurodevelopmental outcomes were computed for patient-related, peri-operative, and post-operative variables and environmental exposure biomarkers. Of the analytes used in this study, several returned nondetectable concentrations for some participants and required imputation consistent with established methods in the literature [29,30]. As such, LOD/sqrt(2) was used to impute values for the two analytes for which the rate of non-detection was low (< 20%) while single imputation was used for the three analytes for which non-detection was more moderate (20 to 50%). For single imputation, specifically, values were drawn from a randomly generated lognormal probability distribution with imputed values bounded between zero and the analyte's LOD. A  total of 453 of 3190 (14.2%) biomarker concentrations were reported as <LOD, and 121 of 3190 (3.8%) concentrations were reported as missing because of insufficient urine volume for analysis or analytical interferences. If the CDC laboratory reported numeric values for concentrations <LOD, these values were used as observed concentrations. To help characterize relationships among the 30 analytes reported here, a table of Spearman-rho correlation coefficients was generated. Because of the extremely large number of bivariate correlations (n = 870) in the matrix, only those achieving statistical significance at an adjusted α = 0.0024 level for 30 moderately correlated endpoints are presented in S2 Table. All observed and imputed urine concentrations were corrected for specific gravity using Levine's formula [31].
Result ADJ ðng=mLÞ ¼ Result � ½ðSG Median À 1:000Þ=ðSG Testing À 1:000Þ� Where Result = observed concentration of an analyte at testing Result ADJ = observed concentration of an analyte at testing after adjusting for specific gravity.
SG Median = median specific gravity value across all subjects. SG Testing = concentration of a subject's urine at the time of testing.

Neurocognitive modeling
Cognitive, Language, and Motor Composite scores from the Bayley Scales of Infant and Toddler Development obtained at 18 months of age were regressed onto a series of 29 patientrelated, peri-operative, and post-operative covariates, individually, using a generalized linear model (normal distribution, identity link) specifically to identify relevant covariates for subsequent modeling of chemical mixtures effects on neurodevelopment. The covariates selected for use in this study were deemed relevant for modeling based upon our prior work, as well as the published literature in CHD [2,32,33]. We elected to keep as many potentially influential covariates as possible for subsequent, in-depth analysis. Assumptions and distributional properties were tested and deemed amenable for analysis. Two of the outcomes were skewed, requiring a Box-Cox transformation to achieve normality and make them amenable for parametric analysis (Cognitive Composite Score^2 . 25 , Motor Composite Score^2 .75 ). Single covariate models with Wald-statistic p-values < 0.20 were used as candidates for inclusion in a best-fitting multiple covariate model for each outcome. Final, best-fitting, multiple-covariate models were then selected based on individual and model-specific likelihood ratio tests, AIC, and BIC values. Covariates used in the final, best-fitting models were then used in subsequent modeling that assessed the neurodevelopmental health effects of the environmental chemical mixture (described in turn).

Modeling the effects of environmental chemical mixtures on neurodevelopment
Because the purpose of this study is to determine if early life chemical exposures are potentially important contributors to our understanding of the mechanisms underlying neurodevelopmental disability in children with CHD, two additional semi-parametric modeling techniques were used to investigate the impact of exposures to environmental chemicals.

Bayesian profile regression (BPR)
BPR is a semi-parametric clustering algorithm that is set in a Bayesian framework with Markov Chain Monte Carlo (MCMC) sampling [34,35]. In this study, BPR was used to identify environmental exposure patterns among the children that could be linked to our study outcomes using conventional linear regression [36]. Advantages of using BPR compared to more conventional clustering algorithms such as K means, Euclidean distance, or nearest neighbor clustering, include a lack of constraint on defining the number of clusters a priori, allowances for categorical covariates on which to base the clusters, its ability to handle and account for correlated, as well as missing data, and a variable selection option to help identify the covariates most responsible for observed patterns of clusters.
Missing exposure (or covariate data) values are accommodated with BPR by denoting each missing value as 'NA' and applying multiple imputation throughout the sampler using the full model information [35,37]. Another particular advantage of BPR is the ability to handle categorical data in the clustering. Hence, if data are highly skewed, they can be reclassified to accommodate the clustering algorithm. In our case, highly skewed data were reclassified into quantiles then fit as categorical covariates, such that those with similarly high levels of exposure (i.e., highest quantile) are clustered together and those with similarly low levels of exposure (i.e., lowest quantile) are clustered together" [34,35].
As mentioned, prior to fitting the model, we first transformed urinary biomarker concentrations into tertiles based on the sample size of our cohort and the highly skewed nature of the biomarker data. Target exposure biomarkers were fit as categorical (discrete) covariates and the variable selection option was utilized to identify exposure biomarkers of interest driving the observed clustering patterns. The PReMiuM package (version 3.2.6) available in the R Statistical Computing Platform (version 4.1.2) was used to implement BPR [34]. Default priors available in the PReMiuM package were used. The outcome was excluded from the cluster analysis so that cluster membership is not influenced by neurodevelopmental outcomes but rather only the co-exposure patterns. In the second step, neurodevelopmental outcomes were analyzed using linear regression models. These outcomes included Cognitive Composite score, Motor Composite score, and Language Composite score. The clusters were fit as categorical variables for adjusted as well as unadjusted linear regression models. In each outcomespecific model, the cluster with the lowest exposure biomarker concentrations was set as a reference group. Clusters with regression coefficient estimates that had 95% confidence intervals that did not overlap with zero were considered significantly associated with the neurodevelopmental outcome of interest.

Weighted quantile sum regression (WQS)
WQS is an alternative statistical method designed to assess associations between correlated joint exposures with a health outcome of interest. With WQS, exposure biomarker concentrations are transformed into quantiles and joined into a weighted index, which both reduces the model parameter dimensions and addresses multi-collinearity. An overall unidirectional effect of the mixture is estimated, with individual exposure biomarkers ranked (or weighted) based on their relative contribution to the weighted quantile index. In our implementation of WQS, we performed 100 bootstrap samples for parameter estimation along with 100 repeated holdouts with cross-validation. Model training and validation were performed using a 40%/60% split of the cohort data. We applied the repeated holdouts because it has recently been shown to improve the stability of WQS estimates in the context of small cohort sample sizes [38]. Since WQS handles missing values by automatically imputing to the lowest exposure group, which would be an inappropriate assumption, only 26 chemical biomarkers were included in the WQS model fitting. We fit WQS using tertiles for ranking of each exposure biomarker and assumed a negative exposure-response for the WQS index since we are interested in adverse effects of the chemical biomarkers mixture on neurodevelopment indicators. A sensitivity analysis was performed which allowed for a positive exposure-response for the exposure mixture; however, there was no evidence of a significant positive effects (data not shown). The health outcome-specific models were modeled as a linear function of the WQS index and covariate adjustment was performed using the same covariates for the linear regression analysis as described above. In our presentation of results, our inference was based on the overall mean WQS index effect estimate along with the corresponding 95% confidence intervals. We evaluated the importance of the estimated weights for each exposure biomarker using a threshold of � 0.038 (1/26 [number of exposure biomarkers]) as this is shown to be a good indicator of variable importance [38,39].

Results
Between September 1, 2011 and August 31, 2015, 140 infants were enrolled in the study. Of these, 110 (79%) returned for the evaluation at 18 months of age (Fig 1). The types of CHD included hypoplastic left heart syndrome or variant (n = 32), transposition of the great arteries (n = 39), and other (n = 39). Patent characteristics and operative management variables are shown in Table 2. The only significant differences were larger placenta weight/birth weight and longer DHCA time in the non-returners.
The mean Cognitive Composite score for the entire cohort (n = 110) was 93.2 ± 12.7. Cognitive Composite scores were more than 1SD below the expected mean for 19 subjects (17%) and more than 2SD below the expected mean for 3 (3%). The mean Language Composite score for the entire cohort was 91.9 ± 18.1. Language Composite scores were more than 1SD below the expected mean for 35 subjects (32%), and more than 2SD below the expected mean for 12 (12%). The mean Motor Composite score for the entire cohort was 92.0 ± 11.7. Motor Composite scores were more than 1SD below the expected mean for 24 subjects (22%), and more than 2SD below the expected mean for 5 (5%). With respect to the multiple covariate models for the Cognitive, Language, and Motor Composite scores, patient factors were more important predictors than operative management variables. (Table 3) Presence of a suspected or confirmed genetic anomaly was an independent predictor of worse performance on all three scores. Older gestational age was a predictor of a higher Cognitive score. Higher SES was associated with better Cognitive and Language Composite scores. Female sex was predictive of better performance for Language and Motor skills. Hispanic ethnicity was associated with worse Motor performance. Among the operative management factors, longer hospital length  of stay for the initial hospitalization was associated with worse performance for the Motor score. No other operative management factors were significantly associated with status at 18 months of age. Urine concentration of exposure biomarkers were measured at 18 months of age; and blood levels were measured for lead, mercury, nicotine, cotinine, and hydroxycotinine. A urine sample was not available for 3 subjects. A blood sample was not available for 20 subjects. Mercury, nicotine, cotinine, and hydroxycotinine were detected in only 1-2 subjects each and  these results were not included in the analyses. A summary of the analyte distributions along with their geometric means, standard errors, and corresponding percentiles is presented in S1 Table. A Spearman-rho correlation matrix of the urinary metabolites concentrations was computed. Of the 870 bivariate correlations, 86 (10%) were highly, statistically significant (p � 0.0024, with many correlations ranging from moderate (r = 0.40 to 0.60) to strong (r > 0.60 to 0.97). Phthalate metabolites in particular were more frequently and strongly correlated with one other urinary metabolite. These observations suggest that there are likely common sources for certain groups of contaminants. In future analyses, we will determine whether specific exposure profiles reflect unique patterns of device use (S2 Table).
Detectable concentrations of some pesticides [4-fluoro-3-phenoxybenzoic acid, trans-3-(2,2-dichlorovinyl)-2,2-dimethylcyclopropane carboxylic acid, and 2-isopropyl-4-methyl-6-hydroxypyrimidine] were found in only 6, 15, and 26 subjects respectively; and were  therefore excluded from analysis. Measurable concentrations of para-nitrophenol, a nonspecific metabolite of parathion were detected in 106 subjects. Parathion is no longer used and there are other substances which can confound the measurements, so para-nitrophenol was excluded from further analysis. Very high concentrations of methyl-and ethyl parabens were detected in some subjects. These chemicals are commonly used in diaper creams. Because of concern over contamination during urine sample collection, these biomarkers were excluded from the analysis. To provide context about the magnitude of exposures compared to the general population, the biomarkers included in the final analysis were compared to available data from two published studies [40,41] (S3 Table). The exposures for the current study were similar to previously reported data.
With BPR, we uniformly analyzed 30 chemical biomarkers, including urinary phenols, parabens, triclocarban and metabolites of pesticides and phthalates, and one blood biomarker (Pb) for the 110 subjects who returned for the 18-month evaluation. (Table 1) Three returnees were removed from the BPR analysis due to lack of any biomarker information, resulting in a final population for profile regression analysis of n = 107. Because profile regression can handle missing values in the exposures of interest, we were able to analyze the data of all 107 participants despite several missing values on a select few biomarkers of exposure (e.g., Pb). Six clusters were identified from BPR analysis of 30 environmental chemical biomarkers. In terms of number of subjects, cluster 1 is the largest cluster (n = 29), followed by cluster 2 (n = 28), cluster 5 (n = 21), cluster 6 (n = 18), cluster 4 (n = 7), and cluster 3 (n = 4). Cluster 1 represents the highest exposure group, with concentrations for 13 phthalate biomarkers that fall within the highest concentration tertile (Fig 2A). Conversely, cluster 2 represents the lowest exposure group with nearly all of the biomarkers involved in the clustering being in the lowest concentration tertile. All other clusters are characterized by elevated concentrations for only a few of the biomarkers involved in the clustering. The variable selection procedure indicated that 15 of the 30 chemical biomarkers drove the observed clustering pattern, 13 of which are phthalate metabolites. The only non-phthalate chemicals included in the clustering were a pyrethroid pesticide metabolite, 3-phenoxybenzoic acid (3-PBA) and a phenol, 2,5-dichlorophenol. There was no statistically significant variability among clusters in terms of study population characteristics (S4 Table). The clear partitioning of phthalate exposure biomarker profiles between the six clusters is shown in the heatmap (Fig 2A), suggesting that the BPR clustering algorithm successfully grouped together individuals with similar exposure biomarker profiles and

Discussion
Exposure to potentially neurotoxic chemicals during susceptible periods of rapid brain growth can have profound effects on neurodevelopment. [4,42] In the current study, we demonstrate that toddlers with CHD are exposed to complex mixtures of chemicals (e.g., phthalates, phenols, parabens, pesticides), some of which are known or suspected EDCs and neurotoxicants, during everyday life. There is significant variability in both the pattern and magnitude of exposures; however, the magnitude of exposure is similar to that reported for infants and toddlers in the general population [40,41]. Greater concentrations of biomarkers of exposure to these chemicals, especially phthalates, are associated with poorer performance for language and motor skills at 18 months of age after adjustment for known risk factors for adverse neurodevelopment outcomes. Many of these chemicals are known EDCs and/or neurotoxicants and we identified a greater adverse effect in girls compared to boys. This study adds to the growing body of evidence of the importance of patient and environmental factors, such as fetal brain development, genetic syndromes, and SES, as determinants of neurodevelopmental outcomes in the CHD population. In order to define factors that may predispose to, or protect against, brain injury and adverse neurodevelopmental outcomes in the CHD population, we must consider the neurodevelopmental exposome, i.e., the totality of exogenous and endogenous exposures from conception onward through adult life [43]. Exposure to phthalates, both in utero and later in life, has been associated with neurobehavioral disability. A systematic review of the literature by Ejaredar and associates found evidence that prenatal exposure to phthalates is associated with adverse cognitive and behavioral outcomes, including lower IQ, and problems with attention, hyperactivity, and poorer social communication [44]. In an important study, Verstraete and colleagues investigated the relationship of circulating phthalate metabolites in a cohort of more than 400 critically ill children (ages 0 to 16 years) to subsequent development of attention deficit hyperactivity disorder (ADHD) compared with normal controls [45]. The phthalate exposure explained half of the risk of attention deficit in the cohort. They concluded that "Iatrogenic exposure to DEHP metabolites during intensive care was independently and robustly associated with the important attention deficit observed in children 4 years after critical illness" [45]. In a study of early life phthalate exposures, Li and coworkers found evidence that childhood exposures to phthalate mixtures was associated with increased behavior problems [46].
In both the cluster analysis and the WQS regression analysis, we found that exposure to mixtures of multiple chemical classes (e.g., phthalates, pesticides, phenols, parabens) was associated with worse outcomes. These classes of chemicals have been associated with adverse neurodevelopmental outcomes in other populations [7,10,47]. The most commonly identified exposure was to phthalates [diesters of phthalic acid] used as solvents and plasticizers (to give flexibility) in many types of products such as food packaging, cosmetics, medical devices, toys, dentures, paints, adhesives, and nail polishes. Phthalates are not covalently bound and thus handling these products can lead to significant exposures [41]. In addition, phthalates and pesticides are commonly found in household dust [11,12]. Prenatal and early childhood exposures to phthalates have been associated with poorer neurobehavioral outcomes [6,48]. After exposure, phthalate diesters are rapidly metabolized and excreted in the urine. Urinary concentrations of these metabolites can be utilized as exposure biomarkers and reflect recent exposures (half-lives of 6 to 24 hours). Di(2-ethylhexyl) phthalate (DEHP) is one of the most commonly used phthalates. Four DEHP metabolites (MECPP, MEHHP, MEOHP, MEHP) were included in the high exposure group (cluster 1); both MEHHP and MEOHP were identified in the WQS regression as important drivers of the adverse exposure-response relationship. Metabolites of several other phthalates were associated with poorer outcomes by both the cluster analysis and WQS regression, including metabolites of di-isobutyl phthalate (MHiBP), benzylbutyl phthalate (MBzP), dibutyl phthalate (MHBP), di-n-octyl phthalate (MCPP), and diethyl phthalate (MEP). Non-phthalate chemicals were also important in the clustering including a pyrethroid pesticide metabolite and 2,5-dichlorophenol. The WQS regression also identified exposure biomarkers of other classes of chemicals (phenols and parabens) as potential drivers of adverse outcomes.
As noted in a recent review, "The staggering majority of epidemiological studies have investigated chemicals in isolation rather than taking into account the totality of chemical exposures, which together may have additive, synergistic, antagonistic, or potentiating effects." [42] A major strength of our study is the use of statistical modeling approaches specifically designed to assess the challenges of analyzing the health effects of multiple correlated environmental exposures. Epidemiologic modeling of the health effects of chemical mixtures, and especially EDCs, must deal with the large number of possible exposures and the potentially highly correlated (multi-collinearity) nature of exposure to multiple EDCs. Large numbers of co-exposures and multi-collinearity make it particularly challenging in terms of being able to disentangle causal effects of particular EDCs within the chemical mixture. Consequently, in this study, we utilized statistical alternatives to multivariable regression that have been developed by others to deal with such challenges [36,39,[49][50][51][52]. Specifically, we applied cluster analysis using BPR as one approach that groups observations together based on their similarity in 'exposure biomarker profiles'. Grouping individuals based on similarity in exposure biomarkers can serve as a type of dimension reduction technique for evaluating health effects of chemical mixtures. Cluster analysis also helps to identify typical exposure biomarker patterns in a study population which can identify sub-populations of individuals who are most vulnerable to multiple exposures and susceptible to adverse health effects. WQS regression is another technique applied in our study that was developed as a way to estimate both an overall effect of the chemical mixture and to derive an exposure-response index weighted by the relative contribution that a particular chemical exposure has to health outcomes of interest (response). This approach therefore helps to identify possible chemicals with worse potential adverse health effects among an array of chemical exposures while also estimating the adverse effects of the mixture. In the current study, we applied BPR (a type of cluster analysis) and WQS to deal with multi-collinearity and for dimension reduction in the exposure-response space. These methods have been applied in other recent studies evaluating the health effects of chemical mixtures in children, including neurodevelopmental outcomes [6, 15-17, 36, 53-55]. In our study, we found these distinct statistical methods to be complementary, insofar as cluster analysis describes sub-populations in terms of their exposures to multiple chemicals and associated adverse neurodevelopmental risks, whereas WQS is concerned with assessing the overall mixture effects on neurodevelopmental risks as well as highlighting specific chemical exposures that may be driving such risks. Another strength is the detailed characterization of the cohort in terms of known and suspected risk factors for poor neurodevelopmental outcomes.
Previous studies have shown that the associations of select chemical exposures with neurodevelopmental outcomes are mixture and sex specific [6,56,57]. The mechanisms resulting in sex differences have not been fully delineated. Many of these chemicals are known EDCs, which can have differential effects on girls and boys due to their different hormone profiles [57]. EDCs are defined as exogenous chemicals that can interfere with any aspect of hormone action [57]. The potential mechanism of EDC mediated neurotoxicity is disruption of endocrine signaling, mimicking the actions of thyroid, glucocorticoid, and gonadal hormones [3,42].
Because of the increasing recognition of the toxicity of EDCs, particularly phthalates, there have been efforts in some countries to use alternative plasticizers, including in medical equipment [58,59]. A variety of alternative plasticizers have been utilized [59]. Studies have shown that there is often lower migration of the alternative plasticizers compared to those manufactured with DEHP, but the toxicity profiles (including impact on neurodevelopment) are not fully known [59,60]. Substitution of alternative plasticizers for phthalates may not always improve outcomes [60]. In a study of NICU patients, use of DEHP-free central venous lines for total parenteral nutrition was associated with a significant decrease in cholestasis [61]. However, in a follow-up to the study evaluating phthalates and ADHD, phasing out DEHP containing devices was shown to reduce DEHP exposure, but there was still significant exposures and the incidence of impaired attention was not reduced [62].
There are limitations to this study. The cohort is relatively small which precluded us from to explicitly testing for statistical interaction effects of specific combinations of exposure biomarkers. While one interpretation of the BPR method could be that interactions are accounted for through effects of latent variables (clustering of exposure profiles), the method is not designed to estimate or test interaction effects between specific chemical exposures of interest. Therefore, BPR at best is an indirect way of assessing the possibility of interaction effects on an outcome. Despite the limitations driven by the relatively small study sample size, leveraging BPR is a strength of this study because it allows us to explore exposure-response relationships given specific combinations of biomarker exposure levels. Another limitation of our application of BPR is the use of "hard" clusters. There is inherent uncertainty in cluster allocation with any type of cluster analysis. While BPR propagates uncertainty by fitting with Markov chain Monte Carlo (MCMC) sampling methods, we focused our epidemiological inference on the deterministic clustering allocation. In other words, we ignored the uncertainty of cluster allocation and instead used the "hard groupings" (or best clustering) as determined with the optimal partition algorithm, "partitioning around medoids" (PAM). Despite this limitation, the hard clustering provides a representative and interpretable cluster allocation that allowed us to model the crude and adjusted associations between exposure profiles and the different neurological outcomes. The assessment of exposure was limited to a single time point and the analyses evaluated cross-sectional exposures and outcomes. While relying on a single exposure biomarker for short-lived compounds certainly leads to exposure misclassification, this misclassification bias will typically be non-differential, which entails effect estimates will be attenuated towards the null. Future studies should include exposure and outcome assessment at multiple time points, including peri-operative, across the life span. Although we assessed multiple classes of chemicals, we did not assess other potentially neurotoxic compounds, including air pollution, flame retardants, and heavy metals other than Pb and mercury. Additionally, our exposure data for Pb included several missing observations. Very little comparative data for the exposure biomarkers are available for 18-month-old children. We also note that two of the clusters (C3 and C4) include small numbers of children. While these two clusters did not show a significant association with neurodevelopmental outcomes, this observation may be driven by the small sample size. Additionally, WQS regression assumes homogeneity for the directionality of the exposure-response relationship. This assumption does not allow for the possibility that some short-lived EDCs could plausibly impart a positive effect on the neurodevelopmental outcomes. Although, our sensitivity analysis suggested there was not a significant positive effect of the mixture. Finally, the neurodevelopment assessment was performed early in life and may not be predictive of later outcomes.
In conclusion, we demonstrate that early life chemical exposures are potentially important factors which should be further investigated as an important contributor to neurodevelopmental disability in children with CHD. In addition, we show that environmental exposure to multiple potentially toxic chemicals and EDCs (e.g., phthalates, pesticides, phenols, parabens) is ubiquitous in toddlers with CHD. There is significant variability in both the pattern and magnitude of exposure biomarkers. Greater concentrations of exposure biomarkers are associated with poorer outcomes for language and motor skills, especially for girls. In order to define factors that may predispose to, or protect against, brain injury and adverse neurodevelopmental outcomes in the CHD population; we must consider not only events that occur during medical care, but the neurodevelopmental exposome, the totality of exogenous and endogenous exposures from conception onward through adult life [43].